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Recent work has provided the means to rigorously determine properties of super-hadronic mat¬ 
ter from experimental data through the application of broad scale modeling of high-energy nuclear 
collisions within a Bayesian framework. These studies have provided unprecedented statistical in¬ 
ferences about the physics underlying nuclear collisions by virtue of simultaneously considering a 
wide range of model parameters and experimental observables. Notably, this approach has been 
used to constrain both the QCD equation of state and the shear viscosity above the quark-hadron 
transition. Although the inferences themselves have a clear meaning, the complex nature of the rela¬ 
tionships between model parameters and observables have remained relatively obscure. We present 
here a novel extension of the standard Bayesian Markov Chain Monte Carlo approach that allows 
for the quantitative determination of how inferences of model parameters are driven by experimen¬ 
tal measurements and their uncertainties. This technique is then applied in the context of heavy 
ion collisions in order to explore previous results in greater depth. The resulting relationships are 
useful for identifying model weaknesses, prioritizing future experimental measurements, and most 
importantly: developing an intuition for the role that different observables play in constraining our 
understanding of the underlying physics. 
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I. INTRODUCTION 

Relativistic nuclear collisions provide a unique opportunity to create and study matter at temperatures and energy 
densities well above the boundary where hadronic degrees of freedom become irrelevant and are replaced with quark 
and gluon degrees of freedom. Whereas the energy density inside a hadron might be on the order of 100 MeV/fm^, 
high-energy nuclear collisions produce mesoscopic regions where the average energy density can exceed 10 GeV/fm^. 
Unfortunately, they are far from ideal environments for extracting the properties of super-hadronic matter. Only the 
asymptotic momenta of outgoing particles are experimentally accessible, so three-dimensional models of the dynamics 
are essential for interpreting results and forming rigorous conclusions about either the nature and properties of 
the matter or about the evolution itself. One serious limitation is that the exact nature of the initial deposition 
and dispersion of energy, momentum, charge, and baryon number in heavy-ion collisions is currently only partially 
understood. Nearly any experimental measurement will be sensitive to these initial conditions as well as to combined 
contributions from the proceeding stages of rapid expansion and cooling through partially equilibrated QGP and 
hadronic phases. This entangles contributions from the various components of the underlying physics and obscures 
the interpretation of many experimental results. 

The development of models that can reliably mimic the various complexities of heavy-ion collisions has been a key 
focus of the field and a necessary step towards condensing experimental results into quantitative statements about 
the underlying physics. Unfortunately, determining physical quantities like the shear and bulk viscosities of QGP is 
non-trivial even with high quality models in place. Experimental data can be largely described using these models but 
the conclusions drawn depend strongly on the validity of the initial state model and other model assumptions. This 
problem is only confounded by the computational complexity of typical heavy-ion collision models which, in many 
cases, makes a full exploration of model parameter space prohibitively expensive and, in turn, encourages simplifying 
assumptions. 

The situation of dealing with very complicated models that map high dimensional parameter spaces to rich heteroge¬ 
nous observational data sets is a common one. It’s faced in disparate fields of science such as biology, econometrics, 
and cosmology. The common solution to this problem is to apply a Bayesian inference approach where a posterior 
distribution over model parameters can be mapped out based on their consistency with a set of observables using a 
Markov Chain Monte Carlo (MCMC) or a related approach [Il|2]. This allows for quantitative conclusions about the 
underlying parameters to be made given the validity of the model and prior assumptions. 

A primary difficulty that has historically limited the applicability of a Bayesian approach in heavy-ion physics has 
been the large number of model evaluations necessary to accurately map out a high dimensional posterior distribution. 
This, coupled with the already computationally expensive nature of the simulations, makes a direct application of 
an MCMC approach essentially impossible with currently available hardware. Recent progress in methodology has 
overcome this difficulty through the use of Gaussian process emulators which can be trained to accurately reproduce 
the results of actual simulated models. The reduced computational load has opened the door to quantitative analyses 
of model parameters which has allowed for systematic constraints of the nuclear equation of state and of the QGP 
shear viscosity las]. 

Stating the single point in the multidimensional parameter space that maximizes the likelihood disregards the 
uncertainty with which the number is stated. In contrast, an MCMC exploration of parameter space provides not 
only the uncertainties for each parameter, but the covariances with other parameters along with the complete shape 
of the posterior likelihood. The methods of [3] indeed provide the complete likelihood distribution, and are now 
extracting the field’s first rigorous quantitative conclusions from heavy-ion collisions at the highest energies. 

Although the Bayesian approach has proved fruitful in nuclear physics, progress in understanding the complexities 
of generated posterior distributions and how they’re driven by the experimental data has been slow. It is relatively 
straightforward to understand how parameters drive observables through varying them and seeing what changes. Go¬ 
ing in the opposite direction, seeing how the measurements drive the inferences about parameters, is more challenging 
for several reasons. For example, changing a certain parameter might readily alter a given observable, but there may 
exist some combination of other parameters that can compensate without affecting the remaining observables. Even 
an observable which does not seem to directly depend on a given parameter, might help constrain other parameters, 
which then helps constrain the first parameter. In deciding which observables to either measure or to better measure, 
one would like to know how the width of the posterior likelihood distribution is affected by reducing the uncertainty 
of a given observable. These relations are critical to gaining insight into understanding not only the degree to which 
model parameters are being constrained, but how and why. Here, we present a new set of techniques for addressing 
all of these questions and then apply the techniques to the 14-parameter analysis first presented in [4] . This analysis 
represents the field’s first quantitative evaluation of much of the field’s consensus-based understanding about which 
observables are truly responsible for addressing the community’s most basic questions about the bulk properties of 
nuclear matter and about the evolution of a high-energy heavy-ion collision. 

In the next two sections we review the model, data and techniques used to generate our previous results in mil], 


3 


presenting a wider range of results from the projections of the MCMC procedure than were presented in [4]. In Sec. IV 


we present techniques for determining linear relations between observables and parameters and between uncertainties 
in observables and widths of the posterior distribution in parameter space. We apply these techniques to the heavy-ion 
analysis described in Sec.sjlljand |ral and focus on the ramifications for extracting the equation of state and viscosity. 
Results are summarized and an outlook is presented in the final section. 


II. MODEL AND DATA OVERVIEW 


The details of the model are expounded in more detail in laa, but we will briefly review the basics and focus 
on describing the 14 parameters varied in this analysis along with the set of observables. Our model consists of the 
generation of an initial state which is then fed into a 2-dimensional hydrodynamics simulation followed by transition 
to a microscopic simulation, known as a hadronic cascade, at a transition temperature of Tq = 165 MeV. An assumed 
symmetry based on an invariance to boosts along the beam axis makes it possible to approximate three-dimensional 
treatment with a two-dimensional model that discretizes the two transverse coordinates. 

Ten model parameters are used to vary the initial state at a time tq = 0.8 fm/c, two are used to determine the 
equation of state in the hydrodynamic stage, and the last two are used to determine the shear viscosity and it’s 
temperature dependence in the hydrodynamic stage. The shear viscosity is parameterized as 


s 



(1) 


where {r]/s)o is the viscosity at Tq and rj' describes the temperature dependence. These two parameters will be 
featured in this paper due to their relatively intuitive meaning and relationship with common observables. Two 
parameters also encapsulate the equation of state. The speed of sound, Cg, as a function of the energy density, e, is 
described as 


Cs(e) =cl{eh) +Q- 
Xo = X' 


Xqx 


^oa; + a;2 + X'2 ’ 
X = ln(e/eft), 


( 2 ) 


where Ch and Cs{eh) are the energy density and speed of sound of a hadron gas at the transition temperature, Tq = 165 
MeV. These quantities are calculated by considering a gas of non-interacting hadrons using the masses and spins of 
particles from the Particle Data Group [5]. All resonances with masses below 2 GeV/c^ were included. With this 
prescription, the equation of state is continuous at Tq. The parameter R describe the behavior of Cs just above Tq 
and the parameter X' provides a scale at which the speed of sound approaches 1/3. Increasing X' lowers the speed 
of sound mainly at high energy density or temperature, and increasing R increases the speed of sound, mainly just 
above Tq. 

Five of the ten initial state parameters apply only to the description of RHIG (Relativistic Heavy Ion Gollider at 
Brookhaven National Laboratory) data. Only data from Au-fAu collisions at lOOA GeV + lOOA GeV is considered. 
The other five described the initial state for Pb+Pb collisions at 1.38A TeV + 1.38A TeV from the LHG (Large 
Hadron Collider at CERN). The initial transverse energy density profile that instantiates the hydrodynamics has the 
form [3] 

~ fwn^wni^') y) “1“ (1 /wn)^sat (^5 : (3) 


where the parameter /^n describes the weighting between two other parameterized forms, ewn which is the wounded 
nucleon form [6], and egat, which describes a form more in line with some ideas of saturation, similarly to [7]. 

ewn{x,y) = (1 - exp(-rB(a:;, y)crsat)), (4) 

^CTgat. 

^Z,^^XXh^TB(x,y) (1 - exp(-T^(x,2/)a,at)), 

^CTsat. 

<^sa.t{x,y) = Z^^-X^^i^}pp^Tmm{x,y) (1 - exp(-Tniax(a:;,2/)o-sat)), (5) 

<^sat. 

Tmin = , T^ax = (Ta + Tb)/2. 

-LA + -Lb 

Here Ta and Tb are the areal densities, the projections of the baryon density of a An or Pb nucleus onto the transverse 
plane, and have dimensions of number per area. In the wounded nucleon model, each nucleon that comes within the 
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Parameter 

Min, Max (RHIC/LHC values) 

Description 

R 

-0.9, 2.0 

Equation of State 

X' 

0.5, 5.0 

Equation of State 

(’?/s)o 

0.02,0.5 

Viscosity at To 

r( 

0.0,3.0 

Temperature dependence of viscosity 

z. 

0.8 / 0.7, 1.25 / 1.4 

Energy normalization ratio 

(7sat(mb) 

22 / 38, 44 / 76 

Saturation cross section 

/wn 

0.0 / 0.0, 1.0 / 1.0 

Weight of wounded-nucleon parameterization 

Fo 

0.2 / 0.2, 1.0 / 1.0 

Initial transverse flow 

/ 

'^XX 

0.0 / 0.0, 1.0 / 1.0 

Initial anisotropy of stress-energy tensor 


TABLE L Fourteen model parameters 


nucleon-nucleon cross section, cr^n = 42 mb at RHIC energies and 73 mb at LHC energies, of one of the other 
nucleons, known as a participant, contributes to the energy density. If the parameter crgat is set equal to the nucleon- 
nucleon cross section then each nucleon can only contribute once and additional collisions do not increase the energy 
density. Relaxing dgat < (^nn allows the particle to contribute multiple times. As cTgat approaches zero the form gives 
binary scaling and the resulting energy density is proportional to T^T^. In the saturated form, the energy density is 
principally determined by the smaller of the two areal densities. Thus, if one nucleon overlaps 5 other nucleons, the 
energy density will be only slightly less than if it overlapped 10. 

The parameter represents the energy per unit rapidity per nucleon collision in a dilute reaction relative to 
the measurement of a pp collision, and both forms satisfy the constraint that for a diffuse overlap of nucleons, i.e. 
Ta^Tb ^ (Tpp the energy density scales as binary collisions, e{x^y) (jpp{dE/dy)ppTATB with 1. However, 

because this is the energy density at tq = 0.8 fm/c, and not the final energy, and because the energy density includes 
longitudinal motion, it differs from the measured transverse energy, and Zg was allowed to vary over a small range 
near unity. For our purposes {dE/dy)^^ was set to 2.69 and 6.0 GeV for the two beam energies, though this number 
is somewhat arbitrary because it is multiplied by Zg. Finally, the fifth parameter describes the initial transverse flow, 
which is assumed to have the form 


rp -^0 r\rp 

loo ^loo 


(6) 


where is the stress-energy tensor. For a traceless stress energy tensor, which would be expected for non-interacting 
gluon fields or collision-less massless particles or for conformal hydrodynamics, and also assuming boost invariance, 
Eo would go to unity [8]. The final parameter adjusted the initial anisotropy of the stress-energy tensor. 


Txx — Tyy — + ‘^^xx)^ — ^(1 '^xx)' 


(7) 


This form is traceless (1/3) 2 3 = P and describes the initial shear, which in the Israel Stewart form of 

hydrodynamics is a dynamical variable which relaxes toward the Navier Stokes value. Parameters are summarized in 
Tabled! 

The model is used to produce a set of 30 observables that correspond to experimental measurements performed 
at RHIC and the LHC. These encompass distilled information describing spectra, elliptic flow, and femtoscopic 
correlations in central and mid-central collisions. These observables were chosen because they all have been shown to 
characterize features of thermalized bulk matter. Fifteen of the observables correspond to Au-j-Au collisions at RHIC 
while the other fifteen correspond to Pb+Pb collisions at the LHC. Observables were taken from two centrality classes, 
the top 0 — 5% centrality cut and the set of collisions in the 20 — 30% centrality cut. It was felt that additional data 
with centrality between these values would be redundant because of the smooth behavior of the observations with 
centrality. Further, it would be difficult to assign uncertainties for observables from more peripheral collisions due to 
the difficulty in justifying hydrodynamic treatments for systems whose size is not larger than a thermal wavelength, 
and where a significant fraction of the transverse profiie is in the corona, which does not fuiiy thermaiize. 

Three ciasses of measurements were considered. The first and most basic is spectra and yieids. These were distiiied 
to four numbers for each centraiity. Three were the mean transverse momenta, pt^ for pions, kaons and protons. The 
average was taken over a finite range of pt , iimited at the iow end by experimentai constraints and cut off at the high 
end to minimize the effects of jets which are non-thermai features outside of this anaiysis. The fourth measurement 
was the muitipiicity of pions within the same pt range. Because chemicai equiiibrium was assumed, rather than 
parameterized, and because baryon annihiiation was not inciuded, the yieids of kaons and protons were not inciuded. 
In [ 3 ] it was shown that modei runs that had the same {pt) had the same spectrai shapes, so no resoiving power was 
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lost by considering only one number to characterize the spectra. Those spectra with the same {pt) as the data also 
provided good descriptions of the experimental spectral shapes US]. 

The second class of observables is comprised of femtoscopic radii extracted from two-pion correlations at small 
relative momentum. Extracting these radii has long been a staple of the field [9]. The gaussian radii are functions 
of momentum, and describe the size and shape of the outgoing phase space cloud of the given momentum. The 
three radii, i?out:^side and i^iong, describe the transverse dimension parallel to the momentum, the transverse size 
perpendicular to the momentum and the size along the beam axis respectively. For this analysis, the radii were 
averaged over the various bins in transverse momentum, so that three observables encapsulated the data for each 
centrality class. 

The final observable was V 2 ^ which characterizes the anisotropic transverse flow which is driven by the elliptical 
shape of the initial transverse profile. The anisotropy, V 2 = (cos 20) in its simplest definition, was only considered in 
the 20-30% centrality bin. This bin was chosen because the model used smooth profiles, generated from the average 
aerial profiles for a given impact parameter, and neglected the lumpy conditions which one would expect from the 
finite number of scatterers in a An or Pb nucleus. Even though this bin is probably the least affected by lumps, the 
experimental value was reduced by a factor of 0.91 to more fairly compare to the model [4]. The next major analysis 
planned by this group will incorporate fluctuating initial conditions, and in addition to more accurately calculating 
V 2 ^ could also address the higher components t’s, r ’4 • • •, which are purely driven by the fluctuations. Because V 2 rises 
nearly linearly with transverse momentum, the value of V 2 was averaged over pt^ with the values weighted by pt^ 
so that bins with higher pt were relatively more important. The linear weighting was motivated by performing a 
principal component analysis of the binned values to find the combination that best captured the variability in the 
model runs. Table [H] summarizes the 30 observables in this analysis. 


III. MARKOV CHAIN MONTE CARLO PROCEDURE AND RESULTS 

The standard approach to determining the likely regions of parameters, x, from comparing model values, 
to experimental values, ^exp, is through MCMC, which provides a sampling of parameters x that are chosen weighted 
proportional to the likelihood. For this study the likelihood is chosen to have a simple Gaussian form, 

r (f ~ ^exp,a) 1 /g\ 

£(a;)~exp|-^--1. (8) 

The uncertainties incorporate both experimental uncertainties and the shortcomings of the model. For instance, if 
the equation of state and viscosity were perfectly described, and if the initial state was parameterized most correctly, 
the model uncertainty describes how accurately one would expect to reproduce the final-state observables given the 
missing physics in the model and the uncertainties in the measurements. In this study, which uses the same model 
output as in [4], the uncertainties aa were all set to 6% of the experimental value. A more accurate determination of 
the uncertainty would require detailed model studies to estimate the impact of missing physics, and a more detailed 
understanding of experimental uncertainties. It would not be surprising to find that a few of these observables are 
understood with slightly better uncertainty or that some are somewhat more uncertain. Such an improvement should 
involve input from both the experimental and modeling communities. In [3] the analysis was repeated with 9% 
uncertainties and the widths of the resulting posterior distributions only increased by ~ 20%. 

Due to the infeasibility of performing millions of full model runs, the MCMC procedure implemented an emulator 
in the place of the full model 013]. The emulator uses an interpolation algorithm to determine the observables from 
1200 full-model runs. The first 1000 parameter sets for the full-model runs uniformly covered the 14-dimensional 
space according to latin hyper-cube sampling. The last 200 parameter sets were chosen to be consistent with the 
likelihood as calculated from the earlier runs and to provide better coverage in the likely region. To improve efficiency, 
the emulator calculates principal components rather than each observable. In this way, those linear combinations of 
observables that stay constant throughout the model runs can be neglected. Thus, rather than evaluating all 30 
observables, only 14 principal components were analyzed for this study. Even that number could be reduced because 
the last several components only varied by a few percent of one value of the experimental uncertainty. 

The MCMC procedure was performed with several million random steps according to the Metropolis algorithm. 
This provides a sampling of the posterior distribution displayed in Fig. In addition to projections onto one 
dimension of the parameter space, the off-diagonal plots show two-dimensional projections. When the elliptic shapes 
of the two-dimensional projections lie at an angle, it shows that a specific linear combination of parameters may be 
well constrained, whereas the orthogonal combination may be poorly constrained. For example, the region of high 
likelihood for the projection onto the plane of two equation of state parameters, one can see that if one increases 
R and X' by similar percentages that the likelihood changes little. Similarly, for the two viscosity parameters one 



6 


observable 

exp. value 

Pt weighting 

centrality 

collaboration 

^2,7r + 7r- 

8.14% 

ave. over 11 pt bins from 160 MeV/c to 1 GeV/c 

20-30% 

STAR [E] 

Rout 

5.28 fm 

ave. over 4 pt bins from 150-500 MeV/c 

0-5% 

STAR [n] 

Rside 

4.81 fm 

ave. over 4 pt bins from 150-500 MeV/c 

0-5% 

STAR [H] 

-Rlong 

5.47 fm 

ave. over 4 pt bins from 150-500 MeV/c 

0-5% 

STAR [TT] 

Rout 

4.27 fm 

ave. over 4 pt bins from 150-500 MeV/c 

20-30% 

STAR [H] 

Rside 

3.99 fm 

ave. over 4 pt bins from 150-500 MeV/c 

20-30% 

STAR [n] 

Rlong 

4.53 fm 

ave. over 4 pt bins from 150-500 MeV/c 

20-30% 

STAR [m 

1 

+ 

494.4 MeV 

0.2 GeV/c<pt < 1.2 GeV/c 

0-5% 

PHENIX [El 

{Pt)K+K- 

796 MeV 

0.4 GeV/c < Pt < 1.6 GeV/c 

0-5% 

PHENIX [12] 

{Pt)pp 

1.135 GeV 

0.6 GeV/c <pt< 2.0 GeV/c 

0-5% 

PHENIX [E] 

1 

+ 

487.5 MeV 

0.2 GeV/c<pt < 1.2 GeV/c 

20-30% 

PHENIX [El 

{Pt)K+K- 

792 MeV 

0.4 GeV/c < Pt < 1.6 GeV/c 

20-30% 

PHENIX [El 

{Pt)pp 

1.111 GeV 

0.6 GeV/c <pt< 2.0 GeV/c 

20-30% 

PHENIX [El 

7T~^7T~ yield 

422 

0.2 GeV/c<pt < 1.2 GeV/c 

0-5% 

PHENIX [El 

yield 

188.7 

0.2 GeV/c <pt <1.2 GeV/c 

20-30% 

PHENIX [E] 

^2,7r + 7r~ 

9.56% 

ave. over 11 pt bins from 0.15 to 1.2 GeV/c 

20-30% 

ALICE [El 

Rout 

5.46 fm 

ave. over 7 pt bins from 200 to 900 MeV/c 

0-5% 

ALICE [H] 

Rside 

5.32 fm 

ave. over 7 pt bins from 200 to 900 MeV/c 

0-5% 

ALICE [E] 

Rlong 

5.72 fm 

ave. over 7 pt bins from 200 to 900 MeV/c 

0-5% 

ALICE [El 

Rout 

4.07 fm 

ave. over 7 pt bins from 200 to 900 MeV/c 

20-30% 

ALICE [E] 

Rside 

4.12 fm 

ave. over 7 pt bins from 200 to 900 MeV/c 

20-30% 

ALICE [E] 

Rlong 

4.41 fm 

ave. over 7 pt bins from 200 to 900 MeV/c 

20-30% 

ALICE [TT] 

+ 

1 

459.1 MeV 

0.1 GeV/c < Pt < 1.2 GeV/c 

0-5% 

ALICE [E] 

{Pt)K+K- 

775 MeV 

0.2 GeV/c < Pt < 1.6 GeV/c 

0-5% 

ALICE [E] 

{Pt)pp 

1.137 GeV 

0.2 GeV/c <pt< 2.0 GeV/c 

0-5% 

ALICE [E] 

1 

+ 

455 MeV 

0.1 GeV/c<pt < 1.2 GeV/c 

20-30% 

ALICE [E] 

{Pt)K+K- 

758 MeV 

0.2 GeV/c < Pt < 1.6 GeV/c 

20-30% 

ALICE [E] 

{Pt)pp 

1.110 MeV 

0.2 GeV/c <pt< 2.0 GeV/c 

20-30% 

ALICE [E] 

yield 

1258 

0.2 GeV/c<pt < 1.2 GeV/c 

0-5% 

ALICE [E] 

7r^7r“ yield 

523 

0.2 GeV/c <pt< 1.2 GeV/c 

20-30% 

ALICE [E] 


TABLE IL Observables used to compare models to data from RHIC (STAR and PHENIX collaborations) and from ALICE 
(at the LHC). To account for non-flow correlations, the value of V 2 was reduced by 9% from the experimental values to account 
for the approximation of smooth initial conditions. 


finds that one can find a good match with {r]/s)o « 0.2 with a modest temperature dependence, r]'. One can also 
match with lower values of {r]/s)o if the viscosity then rises with temperature. Another instance where the posterior 
is off-diagonal is in the projection of {rj/s)^ and /wn- The resulting covariance corroborates the arguments that were 
put forward in [7]. 


IV. SENSITIVITY STUDIES 

The principal goal of this paper is to present various methods for understanding the role certain observables, or 
sets of observables, play in constraining the posterior likelihood of given parameters, or sets of parameters. The 
most straight-forward method to determine the sensitivity is to perform the analysis both with and without a given 
observable, or set of observables. An example of this is illustrated in Fig. The two-dimensional posterior likelihood 
projection for the two viscosity parameters, {r]/s)o and ?]', are redisplayed along with projections where either RHIC 
data or LHC data are ignored. This makes it clear that the LHC data is especially important for constraining the 
temperature dependence of the viscosity, r]'. 

However, the method of repeating the statistical analysis without a given observable can be onerous. If one wishes 
to study the sensitivity to Ny = 30 observables, one would repeat the MCMC procedure Ny times. Here, we present 
methods for finding three simple measures of the sensitivity using the output from a single MCMC trace. That 
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FIG. 1. One- and two-dimensional posterior likelihood projections from the MCMC procedure. The one-dimensional projec¬ 
tions (along the diagonal) illustrate the degree to which an individual parameter is constrained by comparing to data, whereas 
the off-diagonal elements illustrate how some linear combinations of two parameters might be either poorly or well constrained. 
The colored boundaries delineate one-, two- and three-sigma regions, where n-sigma refers to likelihoods of e~^ of the 
maximum likelihood. 


information would be of the form of a list of values, a = 1- • • A^mcmc, where A^mcmc might be of the order of one 
million. For each point one would have the parametes and the observables ^a,a- 

The first measure would simply describe how a given observable, would change due to a change in a parameter 
Xi^ while keeping all other parameters constant. Because is not purely linear, one needs to choose over what 
region the derivative is evaluated. Two choices of interest might be the prior or the posterior. Here, we use single 
brackets, (•••), to denote an average of the prior, and double brackets, ((•••)), to denote averages taken over the 
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FIG. 2. Posterior likelihood distributions for the shear viscosity parameters generated using RHIC and LHC data separately 
(a & b) and together (c). We can see that the viscosity at To is well constrained by RHIC data alone but that there is very 
little constraint of its temperature dependence without the LHC data. This is consistent with the fact that LHC data probe 
higher temperatures. Combining the datasets clearly constrains the data better than a single data set. Panels (d-e) show 
corresponding information for the two equation of state parameters. In this case it is clear that LHC data provide the bulk of 
the resolving power. Again, this is expected because the equation of state is hxed at To = 165 MeV, and LHC collisions probe 
regimes further into the region where the equation of state differs across the parameter space. One can also see that the high 
degree of covariance in both hgures (c) and (f), especially for the covariance between the two equation-of-state parameters in 
(f). This shows that even though neither parameter is particularly well constrained, there is a linear combination that is well 
constrained, while the orthogonal linear combination is poorly constrained. 


posterior. For averages over the prior, one can use the lists from the points at which full-model runs were performed 
for the purpose of training the emulator. If calculations were performed randomly throughout the prior, one could 
consider the covariance {SyaSxi)^ where Sya = ya — (Va) and Sxi = Xi — {xi). One can calculate the partial derivative 
of ya with respect to any parameter by a simple matrix inversion, 

(SyaSxi) = ^ (SxjSxi) (9) 

^ 1 ^) = {SyaSxi){5x5x)^\ 

Here, the brackets around dy/dx note that this is the best slope for this region of parameters space, which in this 
case is the prior. This expression matches the usual least-squares expression for the slope of a line in multidimensions. 
Equation can be easily altered to address how the slope looks when focused on the posterior region by making the 
change (• • •) ^ ((••'))• 


A. Response of extracted parameter values to changes in experimental measurement 


Equation shows how a specific model value, ya^ changes when a specific parameter, is varied. However, 
describing how a small change in an experimental measurement, y^^^\ affects the average value of a parameter, 
{{xi))^ in the average posterior value is a different question. This involves understanding d{{xi))/dy^a^'^\ Here we 
consider a general expression for any posterior-averaged function of x, ((/(x))), which for this specific consideration 
will be f{x) = Xi. 


^ jdxf{x)C{x) 
^ ^y(exp) ^dxC{x) 


( 10 ) 


Both terms on the right-hand-side of Eq. (10) can be calculated from the trace. By setting / = Sxi = (xi — Xi), 
with Xi = {{xi)), the second term vanishes in this case, and if one has an expression for the likelihood, one need 
only average Sxi{l/C)dC over the sampling of points from the MCMC trace to calculate the required result. One can 
invoke the emulator to calculate y{x), and thereby determine both jC and dC. 
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For Gaussian likelihoods Eq. (10) becomes especially simple. In that case 


C 


d 


■A-ij 


i:-J{{5xi5y,)), 



{{SxiSxj)). 


( 11 ) 


Here, 6y can be relative to any point because {{Sx)) = 0, and the derivative in the last line is the average slope for 
the posterior region. Eq. 0 is straightforward to calculate from the output of the MCMC. For the case where E is 
diagonal, 


dxi 


,,(e^P) 


{{5yaSxi)) 


( 12 ) 


To better quantify how a specihc observable constrains a given parameter, one may multiple the expression for 
in Eq. ( |l2[ ) by the measure of how much ya changes throughout the prior, because model values 

of an observable must change within the range of the prior if that observable is to provide resolving power. Figure 
displays this quantity for ah observables. One can also see whether the change is positive or negative, which is not 
obvious and may have an opposite sign compared to dyjdx. 

As examples of the sensitivity analysis one can look at the elliptic flow measured at RHIC and at the LHC. In Fig. 
one can see that changing the measurement of at RHIC strongly affects the extracted value of the viscosity at 
To, {r]ls)o, but has little effect on the extracted temperature dependence, y'. In contrast, the measurement of V 2 at 
the LHC more strongly affects rj'^ while playing a more minor role in determining {rj/s)o. Given the fact that the 
LHC explores higher energy densities, this was expected. One can also see that the constraints on the viscosity were 
driven by measurements of V 2 and the multiplicities, whereas constraints on the two equation-of-state parameters 
were driven by a wider variety of measurements. 


B. Relation between experimental uncertainties and widths of posterior parameter distributions 


We now consider how the uncertainty of a specihed observable affects the width of the posterior parameter distri¬ 
bution, i.e., calculate the derivative 


d 


1 dC 


1 dC 


as 


ab 


Aij = {{5xi5xj--^)) - {{5xi5xj)){{-^^^)) , 


ab 


CdT: 


ab 


(13) 


where the last step followed the steps from the previous subsection used to derive Eq. (10). Equation (13) is fairly 


straight-forward to calculate, but can be a bit cumbersome. To hnd a simpler, though approximate, expression we 
assume that the prior distribution of parameters is Gaussian even though the results shown here ah assumed hard 
cutoffs, and that y{x) is linear, so that the posterior distribution is Gaussian. We also use Gaussian likelihoods for 
comparing to experimental values. When comparing calculations with and without these approximations, results 
changed little. 

With the approximations, the overall likelihood is then of a Gaussian form. 


C ~ exp 




(14) 


5xj = Xj — Xj . 


We have assumed that the parameters Xi have ah been scaled to have the same prior width R. Further, in the 
neighborhood of the maximum likelihood we assume that y behaves linearly. The width of the posterior, can be 
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FIG. 3. Model responses of an observable with respect to a given parameter, {{di/a/dxi)) are displayed in the left-side panel. In 
the right-side panel, the change of the inferred value of a parameter with respect to changes in a measurement, d{{xi)) 
are scaled by the spread of model values throughout the prior, Larger absolute values point to measurements which 

play important roles in constraining that parameter. The signs of the response are not always equal for the corresponding 
derivative in the two plots. 


taken from the Gaussian form, 


f d^x ex.p{—A-j^SxiSxj/2}SxiSxj 
f d^x exY){—AF^^5xidxj/2} 



( 15 ) 
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Our stated goal is to find an expression describing how Aij responds to changes in 


dA 


d 


as 


-1 

ah 


as 




(16) 


ah 


= 2 - 


a 


as 




ah 


t-1 


d^at' 


= -A 


d 




as 


-1 

ah 


^i£ — 


dya 

dXn 


l^\ 

\dxk/ 


Aki- 


(17) 


One can now insert the expression for dy/dx given in Eq. (§ into Eq.s fib\ and [Tt] to find the needed result. When 
working with the hard cutoffs for priors rather than Gaussians, we scale all the parameters to have the same prior 
width, (x^) = in addition to being centered at zero, {x) = 0. If one uses the prior distribution to calculate 
dyjdx^ then the response in Eq. 0 needn’t ever access the experimental value, and one doesn’t need to perform 
the MCMC. 

If one wishes to use the posterior distribution to calculate the derivatives, the form for the response simplifies 


a 


as 


jAij = -{{Sxi5ya)){{5yb5xj)), 


(18) 


ah 


which can also be transformed into an expression in terms of derivatives with respect to E, 


d 


dT. 


-Ain — 


dAii dH 


-1 

cd 


d 


dT^ah dT^ah 


ah 

as 


as, 


ah 


-1 

cd 




-A 


_W-1 W-1 
^ca ^hd ’ 


5E 


ah 


= ^ac {{SycSxi)){{5yd5xj))^^ 


dh 


(19) 


In this analysis the uncertainty matrix is diagonal, T^ah = (^a^ahi and Eq. (19) can be used to calculate the resolving 
power for determining a single parameter Xi due to a single parameter ya. 


Ri\a — 


d 


{{SXiY) d(Ja 
^-1 


A •• 

1 


( 20 ) 


{{5xi 


^-AAyc5xi)){{5yd5xA)T.-,A 


Results for R^a are shown in Fig. 

Even when a specific parameter is not well constrained, a combination of that parameter with another might still 
be well constrained. The two equation-of-state parameters are a good example, and the two viscosity parameters also 
have a strong off-diagonal component to their likelihood projections. These two instances are shown in panels (c) 
and (f) in Eig. If one then wants to understand the contribution of a given observable in constraining the two- 
dimensional projection, one needs to define a quantity which effectively measures the likely area of the two-dimensional 
projection onto the parameters Xi and xj. Here we use the determinant of the two-by-two matrix constructed from 
the ij components of A^j. 

( 21 ) 
( 22 ) 


^ij — A' 


ZZ^JJ 


■ A - ■ A ■■ 

^ZJ Z 1 


Rij;a — ^a ^zj — ^a 
OCFa 


,A^.. 


I ^ ^ A — 2 A ^ A 


Dij would represent the product of the eigenvalues of the two-by-two matrix, or equivalently, the square of the area 
covered by the projection. The way in which Dij changes with respect to a given uncertainty can then be calculated 
with the help of Eq. (19). These sensitivities are illustrated in Eig. This clearly demonstrates that the extracted 
viscosity is strongly affected by the V 2 measurements, and that multiplicities are also important. Other observables 
are of secondary, but not negligible importance. Eor constraining the equation of state, femtoscopic radii seem to 
provide the most resolving power, but all other observables contribute significantly. This underscores the importance 
of simultaneously considering multiple classes of observables to constrain the parameter space. 
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FIG. 4. The resolving power of a particular observable in determining a specific parameter, Ri-a defined in Eq. 
displayed for all observables and parameters. 


(20), is 


V. SUMMARY AND OUTLOOK 

The results here are unprecedented for this field, and for the first time illustrate a systematic method for identifying 
the critical links between parameters and observables. Remarkably, the links identified by the procedure reinforced 
the general knowledge of the field. For instance, it was indeed found that measurements of the elliptic anisotropy 
provide strong constraints of the viscosity. Further, the expectation that RHIC data would play a larger relative role 
in determining the viscosity near Tc, and that the LHC data would play a larger role in determining the temperature 
dependence was confirmed. The relatively large role of the multiplicities was not necessarily expected, nor was the 
fact that other observables provide non-negligible resolving power. 

Whereas there was consensus within the field that V 2 would be the most important observable to determine the 
viscosity, the role of various observables in constraining the equation of state was very much in dispute. Figures 
and all show that the femtoscopic source sizes provide the most resolving power. This validates ideas based on low 
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FIG. 5. The resolving power for determining a pair of parameters, Xj, due to a specific observable ya as defined by Rij-a in 
Eq. (22). The sensitivity for the two viscosity parameters is shown in the upper panel and that for the two equation-of-state 
parameters is displayed in the lower panels. Information about eh viscosity is most strongly determined by measurements of 
V 2 and the multiplicities. The equation of state is most strongly constrained by femtoscopic radii. 


pressures leading to more elongated phase space distributions in the outward direction and that the longitudinal sizes 
would increase if transverse expansion was slower, Additionally, source sizes played a role in measuring 

the final entropy, which is also a measure of the equation of state [18]. The mean transverse momentum had long 
ago been pointed out as being sensitive to the temperature, and therefore the equations of state, m, and even V 2 
had been suggested as a means for extracting early pressure [20]. Therefore, it was not surprising to see the resolving 
power for the equation of state in Fig. [^spread across all observables. 

The techniques presented here could play a pivotal role in determining the direction of future experiments. Before 
embarking on an expensive experimental program to improve the measure of a specific observable, one could check 
to see how that improvement might indeed better determine parameters of greatest interest. For example, one could 
determine whether running an accelerator for a specific projectile target combination, or at a new beam energy, or 
with higher statistics would best provide insight into better determining the equation of state. In many cases it would 
behoove the community to pre-analyze a project with modern statistical techniques before investing the cost and 
manpower for the effort. 
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